module poly_utils_mod

  ! This module provides some utility procedures for dealing with polynomials.

  implicit none

  private

  public calc_monomial
  public calc_deriv_monomial
  public calc_tensor_product_poly_integral_coefs
  public calc_tensor_product_poly_integral_matrix
  public calc_poly_integral_coefs
  public calc_poly_integral_matrix
  public calc_poly_1d
  public calc_poly_deriv_1d

  interface calc_monomial
    module procedure calc_monomial_1d_r8
    module procedure calc_monomial_1d_r16_1
    module procedure calc_monomial_1d_r16_2
    module procedure calc_monomial_2d_r8
    module procedure calc_monomial_2d_r16_1
    module procedure calc_monomial_2d_r16_2
    module procedure calc_monomial_3d_r8
    module procedure calc_monomial_3d_r16_1
    module procedure calc_monomial_3d_r16_2
  end interface calc_monomial

  interface calc_deriv_monomial
    module procedure calc_deriv_monomial_1d_r8
    module procedure calc_deriv_monomial_1d_r16_1
    module procedure calc_deriv_monomial_1d_r16_2
    module procedure calc_deriv_monomial_2d_r8
    module procedure calc_deriv_monomial_2d_r16_1
    module procedure calc_deriv_monomial_2d_r16_2
    module procedure calc_deriv_monomial_3d_r8
    module procedure calc_deriv_monomial_3d_r16_1
    module procedure calc_deriv_monomial_3d_r16_2
  end interface calc_deriv_monomial
  

  interface calc_tensor_product_poly_integral_coefs
    module procedure calc_tensor_product_poly_integral_coefs_1d_r8_1
    module procedure calc_tensor_product_poly_integral_coefs_1d_r8_2
    module procedure calc_tensor_product_poly_integral_coefs_1d_r16
    module procedure calc_tensor_product_poly_integral_coefs_2d_r8_1
    module procedure calc_tensor_product_poly_integral_coefs_2d_r8_2
    module procedure calc_tensor_product_poly_integral_coefs_2d_r16
    module procedure calc_tensor_product_poly_integral_coefs_3d_r8_1
    module procedure calc_tensor_product_poly_integral_coefs_3d_r8_2
    module procedure calc_tensor_product_poly_integral_coefs_3d_r16
  end interface calc_tensor_product_poly_integral_coefs

  interface calc_tensor_product_poly_integral_matrix
    module procedure calc_tensor_product_poly_integral_matrix_1d_r8_1
    module procedure calc_tensor_product_poly_integral_matrix_1d_r8_2
    module procedure calc_tensor_product_poly_integral_matrix_1d_r16
    module procedure calc_tensor_product_poly_integral_matrix_2d_r8_1
    module procedure calc_tensor_product_poly_integral_matrix_2d_r8_2
    module procedure calc_tensor_product_poly_integral_matrix_2d_r16
    module procedure calc_tensor_product_poly_integral_matrix_3d_r8_1
    module procedure calc_tensor_product_poly_integral_matrix_3d_r8_2
    module procedure calc_tensor_product_poly_integral_matrix_3d_r16
  end interface calc_tensor_product_poly_integral_matrix

  interface calc_poly_integral_coefs
    module procedure calc_poly_integral_coefs_r8
    module procedure calc_poly_integral_coefs_r16
  end interface calc_poly_integral_coefs

  interface calc_poly_integral_matrix
    module procedure calc_poly_integral_matrix_r8
    module procedure calc_poly_integral_matrix_r16
  end interface calc_poly_integral_matrix

contains

  subroutine calc_monomial_1d_r8(x, i, val)

    real(8), intent(in ) :: x
    integer, intent(in ) :: i
    real(8), intent(out) :: val

    val = x**i

  end subroutine calc_monomial_1d_r8

  subroutine calc_monomial_1d_r16_1(x, i, val)

    real( 8), intent(in ) :: x
    integer , intent(in ) :: i
    real(16), intent(out) :: val

    val = x**i

  end subroutine calc_monomial_1d_r16_1

  subroutine calc_monomial_1d_r16_2(x, i, val)

    real(16), intent(in ) :: x
    integer , intent(in ) :: i
    real(16), intent(out) :: val

    val = x**i

  end subroutine calc_monomial_1d_r16_2

  subroutine calc_monomial_2d_r8(x, y, i, j, val)

    real(8), intent(in ) :: x
    real(8), intent(in ) :: y
    integer, intent(in ) :: i
    integer, intent(in ) :: j
    real(8), intent(out) :: val

    val = x**i * y**j

  end subroutine calc_monomial_2d_r8

  subroutine calc_monomial_2d_r16_1(x, y, i, j, val)

    real( 8), intent(in ) :: x
    real( 8), intent(in ) :: y
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    real(16), intent(out) :: val

    val = x**i * y**j

  end subroutine calc_monomial_2d_r16_1

  subroutine calc_monomial_2d_r16_2(x, y, i, j, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    real(16), intent(out) :: val

    val = x**i * y**j

  end subroutine calc_monomial_2d_r16_2

  subroutine calc_monomial_3d_r8(x, y, z, i, j, k, val)

    real(8), intent(in ) :: x
    real(8), intent(in ) :: y
    real(8), intent(in ) :: z
    integer, intent(in ) :: i
    integer, intent(in ) :: j
    integer, intent(in ) :: k
    real(8), intent(out) :: val

    val = x**i * y**j * z**k

  end subroutine calc_monomial_3d_r8

  subroutine calc_monomial_3d_r16_1(x, y, z, i, j, k, val)

    real( 8), intent(in ) :: x
    real( 8), intent(in ) :: y
    real( 8), intent(in ) :: z
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: k
    real(16), intent(out) :: val

    val = x**i * y**j * z**k

  end subroutine calc_monomial_3d_r16_1

  subroutine calc_monomial_3d_r16_2(x, y, z, i, j, k, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    real(16), intent(in ) :: z
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: k
    real(16), intent(out) :: val

    val = x**i * y**j * z**k

  end subroutine calc_monomial_3d_r16_2

  subroutine calc_deriv_monomial_1d_r8(x, i, dxn, nx, val)

    real(8), intent(in ) :: x
    integer, intent(in ) :: i
    integer, intent(in ) :: dxn ! x deriv order
    integer, intent(in ) :: nx ! max x order of polynomial
    real(8), intent(out) :: val
    
    real(8) :: cx(-nx:nx)
    integer p
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do

    val = merge( 1._8, merge(0._8,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0)

  end subroutine calc_deriv_monomial_1d_r8

  subroutine calc_deriv_monomial_1d_r16_1(x, i, dxn, nx, val)

    real(16), intent(in ) :: x
    integer , intent(in ) :: i
    integer , intent(in ) :: dxn ! x deriv order
    integer , intent(in ) :: nx ! max x order of polynomial
    real(8) , intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    integer p
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0)

  end subroutine calc_deriv_monomial_1d_r16_1

  subroutine calc_deriv_monomial_1d_r16_2(x, i, dxn, nx, val)

    real(16), intent(in ) :: x
    integer , intent(in ) :: i
    integer , intent(in ) :: dxn ! x deriv order
    integer , intent(in ) :: nx ! max x order of polynomial
    real(16), intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    integer p
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0)

  end subroutine calc_deriv_monomial_1d_r16_2

  subroutine calc_deriv_monomial_2d_r8(x, y, i, j, dxn, dyn, nx, ny, val)

    real(8), intent(in ) :: x
    real(8), intent(in ) :: y
    integer, intent(in ) :: i
    integer, intent(in ) :: j
    integer, intent(in ) :: dxn
    integer, intent(in ) :: dyn
    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    real(8), intent(out) :: val
    
    real(8) :: cx(-nx:nx)
    real(8) :: cy(-ny:ny)
    integer p, q
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do

    val = merge( 1._8, merge(0._8,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._8, merge(0._8,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0)
    
  end subroutine calc_deriv_monomial_2d_r8

  subroutine calc_deriv_monomial_2d_r16_1(x, y, i, j, dxn, dyn, nx, ny, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: dxn
    integer , intent(in ) :: dyn
    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real(8) , intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    real(16) :: cy(-ny:ny)
    integer p, q
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._16, merge(0._16,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0)
    
  end subroutine calc_deriv_monomial_2d_r16_1

  subroutine calc_deriv_monomial_2d_r16_2(x, y, i, j, dxn, dyn, nx, ny, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: dxn
    integer , intent(in ) :: dyn
    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real(16), intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    real(16) :: cy(-ny:ny)
    integer p, q
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._16, merge(0._16,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0)
    
  end subroutine calc_deriv_monomial_2d_r16_2

  subroutine calc_deriv_monomial_3d_r8(x, y, z, i, j, k, dxn, dyn, dzn, nx, ny, nz, val)

    real(8), intent(in ) :: x
    real(8), intent(in ) :: y
    real(8), intent(in ) :: z
    integer, intent(in ) :: i
    integer, intent(in ) :: j
    integer, intent(in ) :: k
    integer, intent(in ) :: dxn
    integer, intent(in ) :: dyn
    integer, intent(in ) :: dzn
    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    integer, intent(in ) :: nz
    real(8), intent(out) :: val
    
    real(8) :: cx(-nx:nx)
    real(8) :: cy(-ny:ny)
    real(8) :: cz(-nz:nz)
    integer p, q, r
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do
    
    cz = 0
    do r = 1,nz
      cz(r) = r
    end do

    val = merge( 1._8, merge(0._8,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._8, merge(0._8,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0) &
        * merge( 1._8, merge(0._8,product(cz(k:k-dzn+1:-1)),k-dzn<0), dzn==0 ) * z**merge(0,k-dzn,k-dzn<0)
    
  end subroutine calc_deriv_monomial_3d_r8

  subroutine calc_deriv_monomial_3d_r16_1(x, y, z, i, j, k, dxn, dyn, dzn, nx, ny, nz, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    real(16), intent(in ) :: z
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: k
    integer , intent(in ) :: dxn
    integer , intent(in ) :: dyn
    integer , intent(in ) :: dzn
    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    integer , intent(in ) :: nz
    real( 8), intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    real(16) :: cy(-ny:ny)
    real(16) :: cz(-nz:nz)
    integer p, q, r
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do
    
    cz = 0
    do r = 1,nz
      cz(r) = r
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._16, merge(0._16,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0) &
        * merge( 1._16, merge(0._16,product(cz(k:k-dzn+1:-1)),k-dzn<0), dzn==0 ) * z**merge(0,k-dzn,k-dzn<0)
    
  end subroutine calc_deriv_monomial_3d_r16_1

  subroutine calc_deriv_monomial_3d_r16_2(x, y, z, i, j, k, dxn, dyn, dzn, nx, ny, nz, val)

    real(16), intent(in ) :: x
    real(16), intent(in ) :: y
    real(16), intent(in ) :: z
    integer , intent(in ) :: i
    integer , intent(in ) :: j
    integer , intent(in ) :: k
    integer , intent(in ) :: dxn
    integer , intent(in ) :: dyn
    integer , intent(in ) :: dzn
    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    integer , intent(in ) :: nz
    real(16), intent(out) :: val
    
    real(16) :: cx(-nx:nx)
    real(16) :: cy(-ny:ny)
    real(16) :: cz(-nz:nz)
    integer p, q, r
    
    cx = 0
    do p = 1,nx
      cx(p) = p
    end do
    
    cy = 0
    do q = 1,ny
      cy(q) = q
    end do
    
    cz = 0
    do r = 1,nz
      cz(r) = r
    end do

    val = merge( 1._16, merge(0._16,product(cx(i:i-dxn+1:-1)),i-dxn<0), dxn==0 ) * x**merge(0,i-dxn,i-dxn<0) &
        * merge( 1._16, merge(0._16,product(cy(j:j-dyn+1:-1)),j-dyn<0), dyn==0 ) * y**merge(0,j-dyn,j-dyn<0) &
        * merge( 1._16, merge(0._16,product(cz(k:k-dzn+1:-1)),k-dzn<0), dzn==0 ) * z**merge(0,k-dzn,k-dzn<0)
    
  end subroutine calc_deriv_monomial_3d_r16_2

  subroutine calc_tensor_product_poly_integral_coefs_1d_r8_1(nx, x0, x1, coefs, mask)

    integer, intent(in ) :: nx
    real(8), intent(in ) :: x0
    real(8), intent(in ) :: x1
    real(8), intent(out) :: coefs(nx)
    integer, intent(in ), optional :: mask(nx)
    
    integer :: msk(nx)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do i = 1, nx
      if (msk(i)==1) then
        count = count + 1
        coefs(count) = (x1**i - x0**i) / i
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_1d_r8_1

  subroutine calc_tensor_product_poly_integral_coefs_1d_r8_2(nx, x0, x1, coefs, mask)

    integer , intent(in ) :: nx
    real( 8), intent(in ) :: x0
    real( 8), intent(in ) :: x1
    real(16), intent(out) :: coefs(nx)
    integer , intent(in ), optional :: mask(nx)
    
    integer :: msk(nx)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do i = 1, nx
      if (msk(i)==1) then
        count = count + 1
        coefs(count) = (x1**i - x0**i) / i
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_1d_r8_2

  subroutine calc_tensor_product_poly_integral_coefs_1d_r16(nx, x0, x1, coefs, mask)

    integer , intent(in ) :: nx
    real(16), intent(in ) :: x0
    real(16), intent(in ) :: x1
    real(16), intent(out) :: coefs(nx)
    integer , intent(in ), optional :: mask(nx)
    
    integer :: msk(nx)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do i = 1, nx
      if (msk(i)==1) then
        count = count + 1
        coefs(count) = (x1**i - x0**i) / i
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_1d_r16

  subroutine calc_tensor_product_poly_integral_coefs_2d_r8_1(nx, ny, x0, x1, y0, y1, coefs, mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    real(8), intent(in ) :: x0
    real(8), intent(in ) :: x1
    real(8), intent(in ) :: y0
    real(8), intent(in ) :: y1
    real(8), intent(out) :: coefs(nx*ny)
    integer, intent(in ), optional :: mask(nx,ny)

    integer msk(nx,ny)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (msk(i,j)==1) then
          count = count + 1
          coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) / (i * j)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_2d_r8_1

  subroutine calc_tensor_product_poly_integral_coefs_2d_r8_2(nx, ny, x0, x1, y0, y1, coefs, mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real( 8), intent(in ) :: x0
    real( 8), intent(in ) :: x1
    real( 8), intent(in ) :: y0
    real( 8), intent(in ) :: y1
    real(16), intent(out) :: coefs(nx*ny)
    integer , intent(in ), optional :: mask(nx,ny)

    integer msk(nx,ny)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (msk(i,j)==1) then
          count = count + 1
          coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) / (i * j)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_2d_r8_2

  subroutine calc_tensor_product_poly_integral_coefs_2d_r16(nx, ny, x0, x1, y0, y1, coefs, mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real(16), intent(in ) :: x0
    real(16), intent(in ) :: x1
    real(16), intent(in ) :: y0
    real(16), intent(in ) :: y1
    real(16), intent(out) :: coefs(nx*ny)
    integer , intent(in ), optional :: mask(nx,ny)

    integer msk(nx,ny)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (msk(i,j)==1) then
          count = count + 1
          coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) / (i * j)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_coefs_2d_r16

  subroutine calc_tensor_product_poly_integral_coefs_3d_r8_1(nx, ny, nz, x0, x1, y0, y1, z0, z1, coefs, mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    integer, intent(in ) :: nz
    real(8), intent(in ) :: x0
    real(8), intent(in ) :: x1
    real(8), intent(in ) :: y0
    real(8), intent(in ) :: y1
    real(8), intent(in ) :: z0
    real(8), intent(in ) :: z1
    real(8), intent(out) :: coefs(nx*ny*nz)
    integer, intent(in ), optional :: mask(nx,ny,nz)

    integer msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (msk(i,j,k)==1) then
            count = count + 1
            coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) * (z1**k - z0**k) / (i * j * k)
          end if
        end do
      end do
    end do
  end subroutine calc_tensor_product_poly_integral_coefs_3d_r8_1

  subroutine calc_tensor_product_poly_integral_coefs_3d_r8_2(nx, ny, nz, x0, x1, y0, y1, z0, z1, coefs, mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    integer, intent(in ) :: nz
    real(8), intent(in ) :: x0
    real(8), intent(in ) :: x1
    real(8), intent(in ) :: y0
    real(8), intent(in ) :: y1
    real(8), intent(in ) :: z0
    real(8), intent(in ) :: z1
    real(16), intent(out) :: coefs(nx*ny*nz)
    integer, intent(in ), optional :: mask(nx,ny,nz)

    integer msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (msk(i,j,k)==1) then
            count = count + 1
            coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) * (z1**k - z0**k) / (i * j * k)
          end if
        end do
      end do
    end do
  end subroutine calc_tensor_product_poly_integral_coefs_3d_r8_2

  subroutine calc_tensor_product_poly_integral_coefs_3d_r16(nx, ny, nz, x0, x1, y0, y1, z0, z1, coefs, mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    integer , intent(in ) :: nz
    real(16), intent(in ) :: x0
    real(16), intent(in ) :: x1
    real(16), intent(in ) :: y0
    real(16), intent(in ) :: y1
    real(16), intent(in ) :: z0
    real(16), intent(in ) :: z1
    real(16), intent(out) :: coefs(nx*ny*nz)
    integer , intent(in ), optional :: mask(nx,ny,nz)

    integer msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(mask))then
      msk = mask
    else
      msk = 1
    endif

    coefs = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (msk(i,j,k)==1) then
            count = count + 1
            coefs(count) = (x1**i - x0**i) * (y1**j - y0**j) * (z1**k - z0**k) / (i * j * k)
          end if
        end do
      end do
    end do
  end subroutine calc_tensor_product_poly_integral_coefs_3d_r16

  subroutine calc_tensor_product_poly_integral_matrix_1d_r8_1(nx, x, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    real( 8), intent(in ) :: x(nx)
    real(16), intent(out) :: A(nx,nx)
    integer , intent(in ), optional :: cell_mask(nx)
    integer , intent(in ), optional :: poly_mask(nx)

    integer cell_msk(nx)
    integer poly_msk(nx)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do i = 1, nx
      if (cell_msk(i)==1) then
        count = count + 1
        call calc_tensor_product_poly_integral_coefs( nx, x(i) - 0.5_8, x(i) + 0.5_8, A(count,:), poly_msk)
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_1d_r8_1

  subroutine calc_tensor_product_poly_integral_matrix_1d_r8_2(nx, x, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    real(8), intent(in ) :: x(nx)
    real(8), intent(out) :: A(nx,nx)
    integer , intent(in ), optional :: cell_mask(nx)
    integer , intent(in ), optional :: poly_mask(nx)

    integer cell_msk(nx)
    integer poly_msk(nx)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do i = 1, nx
      if (cell_msk(i)==1) then
        count = count + 1
        call calc_tensor_product_poly_integral_coefs( nx, x(i) - 0.5_8, x(i) + 0.5_8, A(count,:), poly_msk)
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_1d_r8_2

  subroutine calc_tensor_product_poly_integral_matrix_1d_r16(nx, x, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    real(16), intent(in ) :: x(nx)
    real(16), intent(out) :: A(nx,nx)
    integer , intent(in ), optional :: cell_mask(nx)
    integer , intent(in ), optional :: poly_mask(nx)

    integer cell_msk(nx)
    integer poly_msk(nx)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do i = 1, nx
      if (cell_msk(i)==1) then
        count = count + 1
        call calc_tensor_product_poly_integral_coefs( nx, x(i) - 0.5_16, x(i) + 0.5_16, A(count,:), poly_msk)
      end if
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_1d_r16

  subroutine calc_tensor_product_poly_integral_matrix_2d_r8_1(nx, ny, x, y, A, cell_mask, poly_mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    real(8), intent(in ) :: x(nx)
    real(8), intent(in ) :: y(ny)
    real(8), intent(out) :: A(nx*ny,nx*ny)
    integer, intent(in ), optional :: cell_mask(nx,ny)
    integer, intent(in ), optional :: poly_mask(nx,ny)

    integer cell_msk(nx,ny)
    integer poly_msk(nx,ny)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (cell_msk(i,j)==1) then
          count = count + 1
          call calc_tensor_product_poly_integral_coefs( nx, ny, x(i) - 0.5_8, x(i) + 0.5_8, y(j) - 0.5_8, y(j) + 0.5_8, A(count,:), poly_msk)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_2d_r8_1

  subroutine calc_tensor_product_poly_integral_matrix_2d_r8_2(nx, ny, x, y, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real( 8), intent(in ) :: x(nx)
    real( 8), intent(in ) :: y(ny)
    real(16), intent(out) :: A(nx*ny,nx*ny)
    integer , intent(in ), optional :: cell_mask(nx,ny)
    integer , intent(in ), optional :: poly_mask(nx,ny)

    integer cell_msk(nx,ny)
    integer poly_msk(nx,ny)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (cell_msk(i,j)==1) then
          count = count + 1
          call calc_tensor_product_poly_integral_coefs( nx, ny, x(i) - 0.5_8, x(i) + 0.5_8, y(j) - 0.5_8, y(j) + 0.5_8, A(count,:), poly_msk)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_2d_r8_2

  subroutine calc_tensor_product_poly_integral_matrix_2d_r16(nx, ny, x, y, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    real(16), intent(in ) :: x(nx)
    real(16), intent(in ) :: y(ny)
    real(16), intent(out) :: A(nx*ny,nx*ny)
    integer , intent(in ), optional :: cell_mask(nx,ny)
    integer , intent(in ), optional :: poly_mask(nx,ny)

    integer cell_msk(nx,ny)
    integer poly_msk(nx,ny)
    integer i, j, k, count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do j = 1, ny
      do i = 1, nx
        if (cell_msk(i,j)==1) then
          count = count + 1
          call calc_tensor_product_poly_integral_coefs( nx, ny, x(i) - 0.5_8, x(i) + 0.5_8, y(j) - 0.5_8, y(j) + 0.5_8, A(count,:), poly_msk)
        end if
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_2d_r16

  subroutine calc_tensor_product_poly_integral_matrix_3d_r8_1(nx, ny, nz, x, y, z, A, cell_mask, poly_mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    integer, intent(in ) :: nz
    real(8), intent(in ) :: x(nx)
    real(8), intent(in ) :: y(ny)
    real(8), intent(in ) :: z(nz)
    real(8), intent(out) :: A(nx*ny*nz,nx*ny*nz)
    integer, intent(in ), optional :: cell_mask(nx,ny,nz)
    integer, intent(in ), optional :: poly_mask(nx,ny,nz)

    integer cell_msk(nx,ny,nz)
    integer poly_msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (cell_msk(i,j,k)==1) then
            count = count + 1
            call calc_tensor_product_poly_integral_coefs( nx, ny, nz, &
                                                          x(i) - 0.5_8, x(i) + 0.5_8, &
                                                          y(j) - 0.5_8, y(j) + 0.5_8, &
                                                          z(k) - 0.5_8, z(k) + 0.5_8, &
                                                          A(count,:), poly_msk )
          end if
        end do
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_3d_r8_1

  subroutine calc_tensor_product_poly_integral_matrix_3d_r8_2(nx, ny, nz, x, y, z, A, cell_mask, poly_mask)

    integer, intent(in ) :: nx
    integer, intent(in ) :: ny
    integer, intent(in ) :: nz
    real(8), intent(in ) :: x(nx)
    real(8), intent(in ) :: y(ny)
    real(8), intent(in ) :: z(nz)
    real(16), intent(out) :: A(nx*ny*nz,nx*ny*nz)
    integer, intent(in ), optional :: cell_mask(nx,ny,nz)
    integer, intent(in ), optional :: poly_mask(nx,ny,nz)

    integer cell_msk(nx,ny,nz)
    integer poly_msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (cell_msk(i,j,k)==1) then
            count = count + 1
            call calc_tensor_product_poly_integral_coefs( nx, ny, nz, &
                                                          x(i) - 0.5_8, x(i) + 0.5_8, &
                                                          y(j) - 0.5_8, y(j) + 0.5_8, &
                                                          z(k) - 0.5_8, z(k) + 0.5_8, &
                                                          A(count,:), poly_msk )
          end if
        end do
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_3d_r8_2

  subroutine calc_tensor_product_poly_integral_matrix_3d_r16(nx, ny, nz, x, y, z, A, cell_mask, poly_mask)

    integer , intent(in ) :: nx
    integer , intent(in ) :: ny
    integer , intent(in ) :: nz
    real(16), intent(in ) :: x(nx)
    real(16), intent(in ) :: y(ny)
    real(16), intent(in ) :: z(nz)
    real(16), intent(out) :: A(nx*ny*nz,nx*ny*nz)
    integer , intent(in ), optional :: cell_mask(nx,ny,nz)
    integer , intent(in ), optional :: poly_mask(nx,ny,nz)

    integer cell_msk(nx,ny,nz)
    integer poly_msk(nx,ny,nz)
    integer i, j, k
    integer count

    if(present(cell_mask))then
      cell_msk = cell_mask
    else
      cell_msk = 1
    endif

    if(present(poly_mask))then
      poly_msk = poly_mask
    else
      poly_msk = 1
    endif

    A = 0
    count = 0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (cell_msk(i,j,k)==1) then
            count = count + 1
            call calc_tensor_product_poly_integral_coefs( nx, ny, nz, &
                                                          x(i) - 0.5_16, x(i) + 0.5_16, &
                                                          y(j) - 0.5_16, y(j) + 0.5_16, &
                                                          z(k) - 0.5_16, z(k) + 0.5_16, &
                                                          A(count,:), poly_msk )
          end if
        end do
      end do
    end do

  end subroutine calc_tensor_product_poly_integral_matrix_3d_r16

  subroutine calc_poly_integral_coefs_r8(np, poly, x0, x1, y0, y1, coefs)

    integer , intent(in ) :: np
    integer , intent(in ) :: poly(2,np)
    real( 8), intent(in ) :: x0
    real( 8), intent(in ) :: x1
    real( 8), intent(in ) :: y0
    real( 8), intent(in ) :: y1
    real(16), intent(out) :: coefs(np)

    integer i, j, k

    do k = 1, np
      i = poly(1,k) + 1; j = poly(2,k) + 1
      coefs(k) = (x1**i - x0**i) * (y1**j - y0**j) / (i * j)
    end do

  end subroutine calc_poly_integral_coefs_r8

  subroutine calc_poly_integral_coefs_r16(np, poly, x0, x1, y0, y1, coefs)

    integer , intent(in ) :: np
    integer , intent(in ) :: poly(2,np)
    real(16), intent(in ) :: x0
    real(16), intent(in ) :: x1
    real(16), intent(in ) :: y0
    real(16), intent(in ) :: y1
    real(16), intent(out) :: coefs(np)

    integer i, j, k

    do k = 1, np
      i = poly(1,k) + 1; j = poly(2,k) + 1
      coefs(k) = (x1**i - x0**i) * (y1**j - y0**j) / (i * j)
    end do

  end subroutine calc_poly_integral_coefs_r16

  subroutine calc_poly_integral_matrix_r8(nc, np, poly, x, y, A)

    integer , intent(in ) :: nc         ! Number of cells
    integer , intent(in ) :: np         ! Number of monomials
    integer , intent(in ) :: poly(2,np) ! Monomials degrees (x^?, y^?)
    real( 8), intent(in ) :: x(nc)
    real( 8), intent(in ) :: y(nc)
    real(16), intent(out) :: A(nc,np)

    integer i

    do i = 1, nc
      call calc_poly_integral_coefs(np, poly, x(i) - 0.5_8, x(i) + 0.5_8, y(i) - 0.5_8, y(i) + 0.5_8, A(i,:))
    end do

  end subroutine calc_poly_integral_matrix_r8

  subroutine calc_poly_integral_matrix_r16(nc, np, poly, x, y, A)

    integer , intent(in ) :: nc         ! Number of cells
    integer , intent(in ) :: np         ! Number of monomials
    integer , intent(in ) :: poly(2,np) ! Monomials degrees (x^?, y^?)
    real(16), intent(in ) :: x(nc)
    real(16), intent(in ) :: y(nc)
    real(16), intent(out) :: A(nc,np)

    integer i

    do i = 1, nc
      call calc_poly_integral_coefs(np, poly, x(i) - 0.5_16, x(i) + 0.5_16, y(i) - 0.5_16, y(i) + 0.5_16, A(i,:))
    end do

  end subroutine calc_poly_integral_matrix_r16

  subroutine calc_poly_1d(val,x,nx)
    real   (8) :: x
    integer(4) :: nx
    real   (8) :: val(nx)

    integer :: i
    
    do i = 1, nx
      call calc_monomial(x,i-1,val(i))
    enddo

  end subroutine calc_poly_1d

  subroutine calc_poly_deriv_1d(val,x,dxn,nx)
    real   (8) :: x
    integer(4) :: nx
    integer(4) :: dxn
    real   (8) :: val(nx)

    integer :: i
    
    do i = 1, nx
      call calc_deriv_monomial(x, i-1, dxn, nx, val(i))
    enddo

  end subroutine calc_poly_deriv_1d

end module poly_utils_mod
